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ABSTRACT 

A significant fraction of isolated white dwarfs host magnetic fields in excess of a Mega- 
Gauss. Observations suggest that these fields originate in interacting binary systems 
where the companion is destroyed thus leaving a singular, highly-magnetized white 
dwarf. In post-main-sequence evolution, radial expansion of the parent star may cause 
orbiting companions to become engulfed. During the common envelope phase, as the 
orbital separation rapidly decreases, low-mass companions will tidally disrupt as they 
approach the giant’s core. We hydrodynamically simulate the tidal disruption of plan- 
ets and brown dwarfs, and the subsequent accretion disc formation, in the interior of 
an asymptotic giant branch star. These dynamically formed discs are commensurate 
with previous estimates, suggesting strong magnetic fields may originate from these 
tidal disruption events. 
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white dwarfs 


1 INTRODUCTION 


White dwarfs (WDs) with magnetic field strengths in excess 
of a MegaGauss (MG) are observed in both isolated and bi- 
nary systems (Ferrario et al. 2015). WDs in binaries can ei- 
ther be ’attached” (the companion’s radius extends past the 
equipotential surface and mass transfer is occurring through 
L1) or detached” (the companion’s radius is well within the 
equipotential surface). Between 25% and 36% of WDs in at- 
tached binary systems are highly magnetic (polar or interme- 
diate polar cataclysmic variables) while the remaining sys- 
tems have sub-MG or non-measurable fields (Ferrario et al. 
2015; Pala et al. 2020). In contrast, WDs in detached binary 
systems seldom host strong magnetic fields with only ~10 of 
the over 3000 found to date (Rebassa-Mansergas et al. 2016; 
Parsons et al. 2021). Lastly, about 10% of isolated WDs have 
measured field strengths in excess of a MG, these objects are 
commonly referred to as High Field Magnetic White Dwarfs 
(HFMWDs) (Kawka et al. 2007; Liebert et al. 2003). In each 
of these cases the origin of the WD’s magnetic field strength 
is an area of active research. 

HFMWDs tend to be more massive than their isolated 
non-magnetic counterparts (Liebert 1988). This suggests 
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that HF MWD formation may require some type of mass 
transfer during binary interactions. Additionally, a primor- 
dial field origin in which magnetized white dwarfs are de- 
scendants of Ap/Bp stars is unable to explain the quan- 
tity and field strength distribution of HE MWDs in the solar 
neighborhood (Kawka & Vennes 2004; Kawka 2020). If the 
origin were truly independent of binary interactions, then 
the presence of a detached binary should not affect the inci- 
dence of magnetic fields. However, in contrast to the ~10% 
of isolated white dwarfs that host strong magnetic fields, not 
a single young magnetic white dwarf has been observed in 
a detached binary! (Kawka et al. 2007; Liebert et al. 2003). 
Therefore, the formation channel for strong magnetic fields 
in WDs likely requires binary evolution (Tout et al. 2008; 
Nordhaus et al. 2011). 

Recently, research on binary formation for magnetic 
white dwarfs has focused on pathways that involve common 
envelope evolution (CEE) (Tout et al. 2008; Nordhaus et al. 
2011; Ferrario et al. 2020). In CEE, one of the stellar com- 
ponents experiences radial expansion, creating an envelope 


' The systems in Rebassa-Mansergas et al. (2016) and Parsons 
et al. (2021) all contain older WDs and companions with high 
Roche-lobe filling factors. The latter may indicate a history of 
mass transfer, i.e., not necessarily detached. 
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of loosely bound material. This envelope may directly engulf 
an orbiting companion or cause orbital decay that results in 
a common envelope (CE) phase (Paczynski 1976; Harpaz & 
Soker 1994; Siess & Livio 1999; Nordhaus & Blackman 2006; 
Staff et al. 2016; Belczynski et al. 2016; Wilson & Nord- 
haus 2019, 2020). For a more complete discussion of CEs 
see Ivanova et al. (2013) and references therein. Magnetic 
field generation could occur in the envelope during CEE, or 
during subsequent mass transfer event in the post-CE sys- 
tem, or in a merger if the companion disrupts during the 
CE phase (Tout et al. 2008; Nordhaus et al. 2007, 2011; 
Ohlmann et al. 2016a). 

Double degenerate (WD-WD) mergers are one proposed 
pathway to create isolated HFMWDs (Wickramasinghe & 
Ferrario 2000; García-Berro et al. 2012). In this scenario, 
the resultant object is a massive white dwarf with a hot dif- 
ferentially rotating corona that may amplify magnetic fields 
via a dynamo. This seems unlikely to explain the complete 
distribution of isolated HFMWDs, as the average mass of 
HFMWDs is only about 0.1 Mọ more than the average 
mass of non-magnetic WDs. Furthermore, a small number of 
HFMWDs have been discovered in young stellar clusters. In 
these settings, only the most massive WDs would exist (The 
remnants of the most massive intermediate mass stars), and 
therefore mergers of two of these objects would have resulted 
in Type Ia supernovae (Caiazzo et al. 2020). 

Here, we focus on CE phases that result in the merger 
of a white dwarf with a less massive companion (M-dwarf, 
brown dwarf, or planetary companion) (Nordhaus et al. 
2011; Nordhaus & Spiegel 2013; Guidarelli et al. 2019). Dur- 
ing CEE, the companion does not possess sufficient energy 
to eject the envelope and emerge in a short-period orbit. In- 
stead, the orbital separation decreases until the companion 
is tidally disrupted as it nears the proto-white dwarf core. 
The disrupted material forms an accretion disc that may 
amplify the magnetic field and advect it to the white dwarf 
surface where it could anchor. 

Previously, we showed that such discs survive thermal 
evaporation in the hot stellar ambient medium and remain 
stable for at least ~100 orbits (Guidarelli et al. 2019). How- 
ever, these results assumed a disc formed from disrupted 
material was already present and thus are contingent on 
assumptions about the initial mass and disc geometry we 
chose. To verify that these assumptions are reasonable and 
to obtain further insight, we numerically simulate the dy- 
namics of the tidal disruption events and subsequent disc 
formation. Some work has been done on massive star binary 
mergers with low-mass companions, investigating the long 
term evolution of angular momentum as well as the sur- 
face abundances (Chatzopoulos et al. 2020; Wu et al. 2020). 
Tidal disruption events around white dwarfs have also been 
investigated (Malamud & Perets 2020; Veras et al. 2014). 
However, there have been no simulations investigating disc 
formation and structure from tidal disruption around the 
core of an Asymptotic Giant Branch (AGB) star. 

In this paper we present high-resolution, three- 
dimensional, adaptive mesh refinement (AMR) simulations 
of tidal disruption events involving 10, 20, and 30 My ? com- 
panions inside the core of an AGB star that was evolved 
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from a 2 Mọ zero-age-main-sequence star. We find that an 
accretion disc typically forms within six orbital periods at 
~1IRjo. Approximately, ~60% of the disrupted material or- 
bits in a disk while the rest remains gravitationally bound 
in the tidal tail. Even though the disc dynamically forms, 
its resultant structure is similar to the initial conditions as- 
sumed in our previous work (Guidarelli et al. 2019). 

The structure of the paper is as follows: in Section 2, 
we describe our numerical techniques and initial conditions. 
In Section 3, we analyze the evolution of the global angular 
momentum, present the resultant radially averaged angular 
velocity profile of the disc, quantify the fraction of tidally- 
disrupted mass that forms the disc, analyze the flow prop- 
erties and discuss how our results depend on resolution. We 
discuss future work and conclude in Section 4. 


2 NUMERICAL SETUP 


We perform 3D hydrodynamic simulations of Tidal Dis- 
ruption Events (TDEs) with 10 , 20, and 30 My compan- 
ions around the 0.56 Mo core of an AGB star. In contrast 
to CE simulations that utilize the entire star (Sand et al. 
2020; Chamandy et al. 2019a; Chamandy et al. 2019b, 2020; 
Ohlmann et al. 2016b, 2017; De Marco et al. 2011; Ricker 
& Taam 2008; Ricker & Taam 2012), we focus on the inner 
10!? cm to study the dynamics of the tidal disruption event 
. This results in a higher-effective resolution than what is 
typically achieved in global CE simulations that simulate 
the entire primary star. 

The companions were initialized just outside the tidal 
disruption radius. We define this as the distance of the sec- 
ondary from the center of mass of the primary (AGB core) 
where the companion’s gravitational binding energy is equal 
to the average differential external potential multiplied by 
the diameter of the planet. This separation is given as: 


2(5—n) Meore \ '/? 
Rtidal ~ (> =} Rp , (1) 


where Rp and Mp are the companion radius and mass, 
and Meore is the mass of the proto-WD AGB core. This as- 
sumes that the planet is a polytrope with index n, and the 
companion’s radius is small compared to the separation. 

The companions were given an initial orbital velocity 
that was 65% of the Keplerian velocity at that radius, this 
conserves computational resources as it ensures tidal disrup- 
tion within one orbit. The simulations were run for over 100 
orbits of the inner-most regions of the disc after it formed, 
and at least six orbits of the disc at ~1Rọ (i.e. the approxi- 
mate outer edge of the disc). We investigate the dynamics of 
the tidally disrupted matter, the fraction of initial mass that 
forms the disc, the disc geometry, and the timescale of disc 
formation. Details of the numerical grid, initial conditions, 
modifications to the stellar profile, and other aspects of the 
numerical implementation are described below. 


2.1 Numerical Grid Parameters 


To carry out these simulations we utilize AstroBEAR, an 
Eulerian, three-dimensional (3D) multi-physics code (Cun- 
ningham et al. 2009; Carroll-Nellenback et al. 2013). As- 
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troBEAR employs adaptive mesh refinement (AMR) to dy- 
namically resolve regions of the grid that require higher 
resolution while leaving other regions at lower resolution, 
thereby saving substantial computational costs. 

AstroBEAR solves the following hydrodynamic equa- 
tions: 


2P iY. (py) =0, (2) 

“by +V-(pvv) = —VP — pgc(r)é + pag, (3) 

2E LV (E+P)v] = peel)? +a) y, (4) 

E= ; (v-v)+& (5) 
and 

P=Ak,T, (6) 


where p is the density, t is time, v is the velocity, P is the 
pressure, r is the radial distance to a point particle, ag is 
the acceleration due to self gravity, & is the total internal 
energy density, ñ is the number density of particles, ky is the 
Boltzmann constant and T is the temperature. For these 
simulations we use an ideal gas equation of state with adia- 
batic index, y= 5/3. Lastly, g- is the softened gravitational 
acceleration from a point particle : 


r 2 r>h 


($na) ps 
h2 


h/2<r<h 


(Fy (= +32y)) 
hke 


r<h/2 
(7) 


where h is the softening radius, and y=r/h. Mpp is the mass 
of the point particle which is smaller than Moore. The differ- 
ence in mass Meore — Mpp is taken into account with ag as it 
is resolved on the grid (see Section 2.2.1). 

We employ a base grid of 323 with up to five levels of 
factor two refinement for an effective resolution of 10243. 
This yields a smallest grid cell of ~10? cm, approximately 
the radius of the central white dwarf. The total simulated 
box size is (2 x 10!2cm)?. The spline softening radius of the 
central particle is 10!° cm. We simulate approximately one 
day of physical time, consistent with six orbits at the outer 
radii of the resultant disc and over 100 orbits at the softening 
radius. 


2.2 Initial Condition Generation and Modification 


To generate interior profiles for both the AGB star and the 
companions, we utilize MESA (Modules for Experiments 
in Stellar Astrophysics version 12115) (Paxton et al. 2011, 
2018, 2019). MESA is a spherically symmetric, stellar evolu- 
tion code that is capable of simulating the interior density, 
pressure, temperature and composition profiles of a star over 
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Figure 1. AGB density and pressure profiles before (solid line) and 
after (dotted line) modification. The vertical lines highlight the 
simulation box size, the initial companion separation, the gravita- 
tional smoothing radius, and the smallest grid cell. Note that af- 
ter modification, the smallest grid cell aligns with a much smaller 
density and pressure gradients. 


Figure 2. Density as a function of radius for the planetary com- 
panions, having an effective polytropic index, n~ 1.5. Because the 
profiles are already flat within 10° cm (the smallest grid cell), they 
do not require a point particle or modification. The density cor- 
responding to the flat outer section is the average AGB density in 
the companion vicinity. This section is appended to the profile to 
create an initial co-moving envelope of relatively low density ma- 
terial which reduces the effects of the shock between the orbiting 
companion and stationary ambient material. 


various epochs of its evolution. We calculate the full evolu- 
tion of a primary star with a zero-age-main-sequence mass 
of 2 Meo with solar metallicity. Mass-loss on the Red Giant 
Branch (RGB) is governed via a Reimers mass-loss prescrip- 
tion with a Reimers mass-loss coefficient, ng = 0.7 (Reimers 
1975). On the AGB, we use a Bloecker mass-loss prescription 
with a mass-loss coefficient, Ng, of 0.7 (Bloecker 1995). For 
the simulations reported in this paper, we select the density 
and pressure profiles when the radius attains its maximum 
extent on the AGB. These spherically-symmetric profiles are 
mapped onto the grid and form the initial conditions for the 
AGB ambient. 


P [dynes/cm?] 
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Figure 3. Evolution of the 10 Jupiter mass disc at 6.6, 13.2 and 19.8 hrs. Each panel shows the density distribution in the orbital plane 
indicated by colour. Velocity vectors are shown on all panels with magnitude indicated by colour (in units of the local Keplerian velocity) 
and the projected magnitude by the length. 
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Figure 4. Evolution of the 20 Jupiter mass disc at 6.6, 13.2 and 19.8 hrs. Similar to Figure 3 
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Figure 5. Evolution of the 30 Jupiter mass disc at 6.6, 13.2 and 19.8 hrs. Similar to Figure 3 
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Figure 6. Lower resolution version of our 10M; simulation. Al- 
though more dispersive, the results are morphological similar over 
the simulated time. The arrows represent velocity and use the 
same color scale as in Figures 3-5. 


0 5 10 15 20 


Time (hrs) 


Figure 7. The total angular momentum as a function of time 
within a cylinder of diameter 80% the simulation box side length. 
The ~5% dip in angular momentum at ~6 hr is caused by a small 
amount of material accelerated by the companion that exits the 
cylinder only to slow down and fall back in the the following few 
hours. As indicated by the gray shaded region, At 13 hr, the tidal 
arm passes through the cylinder boundary and causes the total 
angular momentum within the cylinder to decrease accordingly. 


Similarly, we use MESA to generate hydrostatic den- 
sity and pressure profiles for each planet and brown dwarf 
companion. With the built in ’make_planets’ module, one 
can create an initial 1D hydrostatic profile, allow this profile 
to relax and form a core, and expose the profile to expected 
radiation. We initialize our companion with solar metallic- 
ity values of Y = 0.25 and Z = 0.02 and select these profiles 
before cooling and core formation which results in a slightly 
larger estimate for the tidal disruption radius. Note that a 
typical core contributes < 10% of the planet mass. 
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Figure 8. Azimuthally averaged momentum per unit mass perpen- 
dicular to radial position at different times of the 10M; simulation. 
The black curve is the Keplerian velocity profile about a point 
particle of mass 0.565 Mo (Mpp). While our simulations result in 
elliptical discs, we expect circularization in a few more orbital 
time scales. When this occurs, the disc will still be sub-Keplerian 
as pressure gradients will offset some gravitational force. 


2.2.1 AGB Profile Modification 


Because of the large range of scales inherent to AGB stars, 
we utilize a point particle to represent some of the mass 
in the central region. Were a point particle not employed, 
the steep density gradients in this region would be under- 
resolved, resulting in numerical instability on small scales 
and an inaccurate gravitational potential at larger scales. 
Once implemented, point particles allow for the original 
steep density profile to be replaced with a flatter profile. 
However, point particles also require a modification to the 
gravitational potential as under-resolved small separations 
lead to extraneous acceleration. Therefore, to maintain a 
proper representation of the potential when implementing a 
point particle, we modify the density profile in accordance 
with a smoothed gravitational potential acting on the sur- 
rounding material. This modified profile is shown in Figure 
1. 

We modify the density and pressure distributions using 
a method previously described in Guidarelli et al. (2019). 
Specifically, we construct a numerical integrator to iterate 
over solutions to the Modified-Lane-Emden equation. This 
yields a modified density and corresponding pressure pro- 
file that maintains hydrostatic equilibrium for any choice of 
initial profile and smoothing function. This code is publicly 
available? and is useful for generating hydrostatic initial con- 
ditions for 3D simulations involving evolved stars (Guidarelli 
et al. 2019; Ohlmann et al. 2017; Chamandy et al. 2018). 


2.2.2 Planet and brown dwarf co-moving envelopes 


In contrast to the AGB initial conditions, the planetary 
profiles could be imported to the 3D grid without mod- 
ification, as the smallest grid-cell accurately resolves all 
density gradients. However, because the pressure gradients 
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Figure 9. Contours of constant density (0.02 g/cm?) in planes perpendicular to the orbital plane at 20 hrs into the 10 My simulation. Color 
indicates the azimuthal angle in increments of 45 degrees about the point particle. Superimposed is a to-scale image of the 10 My disc 
from our previous work (Guidarelli et al. 2019). The red blob on the left is a cross-section of the tidal tail. 
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Figure 10. Vertical slice through the orbital plane showing the 
internal temperature of the disc at 26.4 hrs into the 10 My com- 
panion simulation. The disc is much colder than the stellar inte- 
rior and entrainment of the of ambient material is caused by the 
intense vertical shear. 


at the planetary boundary cause large initial shocks, we 
give the surrounding ambient (a spherical shell of thickness 
0.5 Rp) the same initial velocity as shown in Figure 2. While 
shocks are expected, the initial steep pressure gradient sig- 
nificantly lengths the computational time. By adding this 
low-density shell of co-moving material, the initial shock dis- 
sipates quickly allowing subsequent shocks to develop natu- 
rally. 


3 RESULTS 


Figures 3 - 5 present the density distributions in the orbital 
plane from each simulation at 6.6, 13.2 and 19.8 hrs. These 
times correspond to ~2, 4 and 6 Keplerian orbits at the 
outer disc radii. In each panel we also show the velocity in 
units of the local Keplerian circular speed. From these im- 
ages, it is clear a disc begins to form within 20 hrs regardless 


of the initial companion mass. Within the first few hours the 
planets disrupt, dispersing the angular momentum. The ma- 
terial that gains angular momentum forms a tidal tail and 
moves outward while material that loses significant angular 
momentum falls inward. The main differences between these 
simulations stem from the increase in total linear momen- 
tum. During disruption, the more massive companions trans- 
fer more linear momentum to the primary core. Because the 
core then moves in the direction of this momentum, material 
in the opposite direction becomes more dispersed. 

Thermal evaporation seems to play a minor role as the 
disrupted companion mass remains relatively distinct from 
the hot stellar interior. At the end of the simulations, ~60% 
of the initial companion mass resides in the disc. The ma- 
terial in the tidal tail remains bound and may fall back on 
longer timescales. 

To quantify the effect of resolution on our results, we 
present a snapshot from an identical 10 My simulation at 
lower resolution in Figure 6. Although more dispersive, the 
lower resolution simulation is dynamically and morphologi- 
cally consistent with the higher resolution run. 


3.1 Angular Momentum Evolution 


The total angular momentum within a cylinder of radius 
0.8 x 10!2cm is shown in Figure 7 as a function of time for 
each of the simulations. In all cases angular momentum is 
conserved to within ~5% until material from the disrupted 
companion leaves the computational grid at ~13hrs. 

The distribution of angular velocity for the 10 My planet 
is shown in Figure 8. Each color represents a different time in 
the simulation in reference to the Keplerian velocity (black 
line). As the planet moves inwards, density and angular mo- 
mentum diffuse radially. By the end of the simulation, the 
angular velocity profile matches the concavity of the Kep- 
lerian solution. The difference in magnitude is the result of 
a non-circular (elliptical) disc as well as pressure support 
within the disc. 
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3.2 Resultant Disc Geometry and Efficiency 


Figure 9 presents density contours of the disc taken from 
slices through the point particle perpendicular to the orbital 
plane at various angles (indicated by color). In this figure, 
we included an image of a 10M; disc from our previous work 
for comparison (Guidarelli et al. 2019). In aggregate, the 
contours provide an estimate of the height to radius ratio of 
the disc. We find that for these discs H/R varies from 0.05 
in the inner regions to 0.25 in the tidal tail consistent with 
our previous work (Guidarelli et al. 2019). The temperature 
in a vertical slice through the disc is shown in Figure 10 and 
demonstrates that post-tidal disruption, the disk remains 
distinct and cold compared to the hot ambient. 


4 CONCLUSIONS 


In this work, we present the results of 3D hydrodynamic sim- 
ulations of tidal disruption events of planet and brown dwarf 
companions in the interior of an AGB star. These events are 
expected to form accretion discs which may result in condi- 
tions that amplify and anchor strong magnetic fields to the 
proto-WD (Nordhaus et al. 2010). As mass-loss from the 
AGB star continues during its evolution, a HFMWD would 
eventually emerge. We show that for planetary and brown 
dwarf companions, tidal disruption results in the formation 
of an accretion disc around the proto-WD which typically 
occurs within a few orbital timescales. About 60% of the 
original companion mass forms the disc with the remain- 
ing 40% in the tidal tail and subject to fallback on longer 
timescales. Lastly, the innermost region of our tidal discs are 
morphologically similar to that of the disc initial conditions 
assumed in our previous work (Guidarelli et al. 2019). 

Once the disc is formed, advection of magnetic flux 
through the disc or its surface might transport the field to 
the white dwarf. Note that this would occur on radial accre- 
tion timescales* (< 1 month) for a standard disk, which is 
much shorter then a typical AGB lifetime (~10°~® years). 

Within the environment of a stellar envelope however, 
a buildup of pressure may halt accretion in the absence of 
a pressure release valve. A pressure release valve could be 
provided by jets (Chamandy et al. 2018). The magnetic field 
itself may incite vertical outflows which act to relieve pres- 
sure, in essence forming a self-sustaining conveyor belt of 
magnetic flux until the field gets large enough to inhibit 
further advection. Future work should therefore include a 
study of the dynamics of tidal disruption when the compan- 
ion, and/or stellar ambient, are initially magnetized. 
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